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Abstract 

In this paper the traveUing wave solutions in the adiabatic model with 
two-step chain branching reaction mechanism are investigated both nu- 
merically and analytically in the limit of equal diffusivity of reactant, 
radicals and heat. The properties of these solutions and their stability are 
investigated in detail. The behaviour of combustion waves are demon- 
strated to have similarities with the properties of nonadiabatic one-step 
combustion waves in that there is a residual amount of fuel left behind the 
travelling waves and the solutions can exhibit extinction. The difference 
between the nonadiabatic one-step and adiabatic two-step models is found 
in the behaviour of the combustion waves near the extinction condition. It 
is shown that the flame velocity drops down to zero and a standing com- 
bustion wave is formed as the extinction condition is reached. Prospects 
of further work are also discussed. 

Keywords: combustion waves, Evans function, chain-branching reaction, 
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1 Introduction 

Combustion waves have been studied for some time and are topic of a relatively 
recent review [I]. They have been observed in numerous experiments T and 
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play an important role in industrial processes, such as the production of ad- 
vanced materials using Self-propagating High-temperature Synthesis (SHS) [2]. 
The one-step irreversible reaction models have contributed greatly to our under- 
standing of these phenomena. In these models it is assumed that the reaction is 
well modeled by a single step of fuel (F) and oxidant (O2) combining to produce 
products (P) and heat. The generic kinetic schemes of models with one-step 
chemistry are: F — > P + heat oi F + O2 — > 2P + heat, where the temperature 
dependent rate of the reaction is given by Arrhenius kinetics k(T) = e""^"/"^ 
and Ta is the activation temperature. These models have proven their useful- 
ness since they are relatively simple and allow analytical investigation using 
asymptotic methods in the limit of infinitely large activation temperature [31 [3] . 
The most important feature of one-step models is that they have led to many 
useful and qualitatively correct predictions for phenomena such as: ignition, ex- 
tinction and stability of diffusion flames; propagation and stability of premixed 
flames; flame balls and their stability; structure and propagation of flame edges 
etc. 

However, in the overwhelming majority of cases, the chemical reactions in 
flames proceed according to a complex mechanism, that involves a variety of 
different steps [5]. Moreover, for many reactions, models with simple one-step 
kinetics may lead to erroneous conclusions as noted in [S]. In other words, if 
we want to obtain a realistic description of the flame kinetics, several differ- 
ent steps, each with its own intermediate chemical species, have to be taken 
into account. Recent advances in computational power have made it possible 
to study the flame behaviour using full numerical solutions of the equations of 
energy and mass transfer for all of the species involved with detailed chemistry. 
Several numerical codes have been developed in order to carry out these calcula- 
tions, such as Sandias CHEMKIN code [B] or the Flame Master code [7j. These 
numerical algorithms have been successfully applied to analyze the properties 
of both diffusion [8j and premixed flames [9] . Although such investigations are 
useful in providing some quantitative prediction for observed phenomena, there 
is still a great deal of uncertainty about the reliability of these complex models 
when applied to the prediction of generic behaviour of flames such as stable 
combustion regimes, limits of the flame extinction and particularly the onset 
of combustion phenomena such as pulsating and cellular combustion, when the 
reactions begin to change rapidly in space and/or time. 

A logical compromise between the models with singe-step and detailed multi- 
step kinetics have been found recently in using reduced kinetic mechanisms. In 
a number of papers [51 [HI [THl dH (listing only few of them) the detailed schemes 
of the hydrogen and methane oxidation, which include dozens of intermediate 
reactions, were successfully reduced to several steps. The remarkable feature 
of these type of models is that they allow on one hand, analytical investigation 
[TUl [TT] to be successfully undertaken, while on the other hand they are able to 
produce excellent quantitative results jSKS] to accurately predict the main flame 
characteristics for some specific reactions such as flame propagation velocities, 
flame structures including the profiles of temperature and reactants etc. These 
studies are very important for applied problems, in which the properties of 
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flames for specific reaction and flame conflguration are studied. 

The importance of investigation of flames with complex kinetics is wefl ac- 
knowledged in combustion literature and as noted in Clavin et al. [T2] "a crucial 
first step in such a study consists in determining reduced schemes that capture 
the essence of complete network" . A particulary promising candidate is a two- 
step reaction mechanism. In this paper we only consider combustion waves in 
premixed fiames. For premixed fiames, several models involving two chemical 
steps have been considered previously. These studies address the fundamen- 
tal (generic) properties of fiames with multi-step reaction mechanisms rather 
than consider some specific combustion reaction. The two-step reactions can be 
split into two groups: reactions with competing steps [31 [T^ and sequential steps 
[3l[T3l[T4l[T5l[T6l[T7l[T8j. In the first group of models three types of reactions are 
considered: (i) simple competing scheme Ci: A Pi, A —* P2; (ii) competing 
chain reaction scheme C2: A ^ B, A + B P; (iii) competing chain branching 
scheme C3: A + B ^ 2B, A + B ^ P, where A is the reactant, B the radical, 
P the product. The second group of the reactions includes several schemes: 
(i) non competitive reactions Ni: [TB] Ai ^ Pi, A2 ^ P2', (ii) sequential fuel 
decomposition reaction N2: [H [Ml [17] A ^ B ^ P; (iii) chain branching re- 
action N3: ^\Wl\M\IB\M ^ + B -^2B,nB + M ^nP + M, where n = 1 
in [131 [TS] and n = 2 in [31 [13] and M is a third body. These schemes were 
investigated either analytically by using high activation energy asymptotic or 
numerically. As a rule the asymptotic methods allow the investigation only in 
certain distinguished limits when the inner equations arising in the asymptotic 
procedure can be integrated analytically [12l [HI [15l [17]. However, generally 
these equations cannot be integrated analytically due to the increased complex- 
ity of the two-step models in comparison to one-step reaction models and the 
inner equations are solved numerically to match them to the outer equation 
solutions [ill [HI [13] ■ In most cases the two-step reaction models show more 
rich dynamical behaviour in comparison to their single step counterparts. For 
given parameter values there can exist combustion wave solution branches trav- 
elling with unique speed, with two different speeds (C-shaped dependence of 
speed on parameters), and with three different speeds (S-shaped dependence 
of speed on parameters). The multiplicity of the flames suggests that a more 
complex stability behaviour also exist. However, the stability analysis of these 
flames are reported in only a few papers |121 115j and has not been investigated 
systematically. Since the behaviour of multi-step models can differ from their 
one-step prototypes, we can expect a number of new phenomena such as bista- 
bility (which cannot be predicted by the one-step models) to be found as a result 
of such an analysis [12]. 

In this paper we focus on the investigation of premixed combustion waves 
in a model with non-competing two-step chain branching reaction. As noted 
in |14| the chain branching reaction is more realistic in describing real flames 
such as hydrocarbon flames in comparison to other two-step reaction models. 
This model was initially considered in [3J and then generalized in |13| . It is 
usually referred to as Zeldovich-Liiian model. In [T3] the model is considered 
in a limiting case when the recombination step is fast in comparison to heat 
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conduction time. Firstly,the resulting problem is investigated semi-analytically 
by using activation energy asymptotic first which provides a set of equations 
for the inner problem. These equations were then solved numerically. It was 
shown that there exists a unique travelling wave solution for given parameter 
values. In [TT] a reduced two-step chain branching reaction for oxidation of 
hydrogen was studied. The methodology used in this paper is similar to [T3j (i.e. 
semi-analytical investigation using activation energy asymptotics and numerical 
integration of the resulting equations). It was shown that the flame has a unique 
burning velocity. The results obtained in this paper are shown to be in a good 
agreement with numerical experiments for the detailed kinetic mechanism of 
the reaction, which indicates that the two-step reaction models might have 
a capacity of quantitatively predicting the flame behaviour. There are other 
papers devoted to investigation of the premixed flames in this model and we 
refer the reader to |15j for more detailed overview, however the stability of 
combustion waves in the Zeldovich-Lihan model has not been investigated. In 
[m [13] a simplified version of Zeldovich-Lihan model is introduced. In this 
model the order of the recombination reaction is taken as n = 1. In contrast 
to [3] the model also includes the linear heat loss to the surroundings. The 
model is studied in the limit of high activation energy which allows asymptotic 
analysis to be carried out successfully. The speed of the combustion wave is 
determined as a function of the parameters of the problem. It is shown that in 
the nonadiabatic case the dependence is C-shaped. For the adiabatic case the 
expression derived in [15] suggests a unique speed of the flame propagation. It is 
remarkable that in [15j the stability analysis of combustion waves is carried out 
for the case of the reactant Lewis number less than one. The analysis predicts 
that the wave can lose stability due to the cellular instabilities emerged in this 
case. In our previous paper [^Oj we have investigated the properties of the model 
introduced in [15j in the adiabatic case and in the limit of equal diffusivity of 
the reactant, the radical and heat. In contrast to [15] the activation energy in 
[20] is taken to be an arbitrary number (not an infinite number) and we used a 
different nondimcnsionalization. Our non dimensional parameters are common 
for one-step models, which makes for an easier comparison between the two 
models. In [50] the properties of the travelling wave solutions are investigated 
in detail by means of numerical simulation. It is demonstrated that the speed 
of combustion wave as a function of parameters is single valued. We have found 
that for finite activation energy there is a residual amount of reactant left behind 
the travelling combustion wave which is not used in the reaction. This makes 
the problem similar to the nonadiabatic one-step premixed flames. The other 
fact about the model considered in [20j which makes the similarity between the 
adiabatic two-step reaction and the nonadiabatic one-step system even stronger 
is that, at certain parameter values the combustion wave exhibits the extinction. 
However, it occurs in such manner that the speed of the wave drops to zero. 
This mainly distinguishes the one- and two-step models, since in the former case 
the wave can propagate for any parameter values and even in the nonadiabatic 
case the speed of combustion wave does not vanish for finite parameter values, 
whereas this is not the case for the two-step models. 
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Extinction of the combustion wave is a critical property of the model that 
also distinguished the two-step models form their one-step counterparts. There- 
fore we have also examined this property in detail in this paper. The other 
important question which will be studied here is the stability of the combustion 
wave. We will compare our results (where possible) with the predictions of the 
asymptotic analysis of (TS]. In [50] the existence of multiple travelling wave 
train solutions is reported for some parameter values. Our analysis shows that 
besides wave trains, the solutions with various structure can also exist such as 
pulses, fronts with step-like structure etc. The properties and relevance of these 
solutions lies beyond the scope of this paper and will be discussed in future 
work. In this paper we will only consider the classical travelling front solutions. 

The rest of the paper is organized as follows. The nondimensional model 
equations are introduced in section [2l and where possible, we have referred the 
reader to [20j to avoid repetition. In this section we also outline the formula- 
tion of the problem to obtain the travelling front solution and determine the 
condition for the existence of the travelling waves. Section 3 is devoted to the 
investigation of the travelling wave solution. In the first part of section 3, we 
summarize the properties of this solution which were obtained in [20] . We then 
determine the extinction limit for the travelling front solution. We analyze its 
behavior near the point of extinction using the matched asymptotic expansion 
method and compare the predictions of analytical and numerical approaches. 
We conclude section 3 by undertaking the stability analysis of the travelling 
wave solutions. Finally, in section 4 the concluding remarks are presented with 
some discussion for future work. 



We consider an adiabatic model (in one spatial dimension) that includes two 
steps: autocatalytic chain branching A + B ^ 2B and recombination B + M ^ 
C + M. Here A is the fuel, B the radical, C the product, and M a third body. 
Following [31[Tll|Tn] we assume that all the heat of the reaction is released during 
the recombination stage and the chain branching stage does not produce or 
consume any heat. As noted in [3], in this scheme the recombination stage serves 
both as an inhibitor which terminates the chain branching and an accelerant 
which produces heat. According to [20]) the equations governing this process 
can be written in nondimensional form as 



where u, v and w are the non-dimensional temperature, concentration of fuel 
and radicals respectively; x and t are the dimcnsionless spatial coordinate and 
time respectively. In ([T]) the following non-dimensional parameters have been 
introduced: f3 — CpE/QR, r ~ Ar/As, ta.b = pCpDA.s/k, where Cp is the spe- 
cific heat; Da and Db represent the diffusivity of fuel and radicals respectively. 



2 Model 



Ut = Uxx + rw, 
vt = TAVxx - (3vwe~ 
wt = tbWxx + Pvwe 




r/3w, 



(1) 
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Ar and A b are constants of recombination and chain branching reactions respec- 
tively; Q is the heat of the recombination reaction; E is the activation energy for 
chain branching reaction; R is the universal gas constant. Here parameter (3 is 
the dimensionless activation energy of the chain-branching step which coincides 
with the corresponding definition for the one step model [3] . Parameter r is the 
ratio of the characteristic time of the recombination and branching steps and 
cannot be reproduced in one step approximation of the flame kinetics. 
Equations ((T|) are considered subject to the boundary conditions 

u = 0, V = I, w = for X — > oo, , , 

Ux ~ 0, Vx — 0, w ^ for X — > — oo. 

On the right boundary we have cold (u = 0) and unburned state (v — 1), where 
the fuel has not been consumed yet and no radicals have been produced {w = 0). 
We also take the ambient temperature to be equal to zero. As noted in [19] this 
is a convenient way to circumvent the so-called "cold-boundary problem" and 
it does not change the generic behaviour of the system. On the left boundary 
(x — > — oo) neither the temperature of the mixture nor the concentration of fuel 
can be specified. We only require that there is no reaction occurring so the 
solution reaches a stationary point of ([T]). Therefore the derivatives of u, v are 
zeros and w — for x — > — oo. 



2.1 Travelling wave solution 

We seek a solution to the problem ([1]), ((21) in the form of a travelling wave 
u{x,t) — u(^), v{x,t) — f(^), and w{x,t) — w{S.), where we have introduced 
^ = a; — ci, a coordinate in the moving frame and c, the speed of the travelling 
wave. Substituting the travelling wave solution into the governing equations we 
obtain 

Mjj -|- cu{ + rw = 0, 

TAV^^ + cv^ - /Jwwe"!/" 0, (3) 
tbw^^ + cw(^ + Pvwe~^^" — r[3w = 0, 

and boundary conditions 

u = 0, V = 1, w = Q for ^ oo, , 
= 0, 0, w = for ^ — oo. ^ 

Following [21 [inj we consider the case when Lewis numbers for the fuel and the 
radicals are equal to unity. Although this assumption simplifies the problem 
significantly, it still allows the investigation of the generic properties of the 
system © and (g]). 

In the case ta = tb = ^ equations ^ possess an integral S = (3u-\- v + w. 
Using S and the first boundary condition in ([4]), equations Q can be reduced to 
a system of two second order equations for temperature and fuel concentration 

-I- cu^ + r(l — /3u — = 0, , , 
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Figure 1: Schematic representation of the traveUing wave solution on (/3m, v) 
plane. 

On the right boundary we require that u = and v = 1, whereas on the left 
boundary — > — oo) we modify the boundary conditions as follows 



where a denotes the residual amount of fuel left behind the wave and is unknown 
until a solution is obtained. We note here that at first glance, system ([5]) looks 
very similar to the equations describing the dynamics of the one-step adiabatic 
reaction model considered in [19j . However, in contrast to the one-step adiabatic 
case, equations ([5]) do not have an integral which enabled further simplification 
in [12]. Moreover, boundary conditions (O suggest that there can be some 
fuel left behind the reaction zone, which is impossible in the case of a one-step 
adiabatic reaction model. 

2.2 Existence of the travelling wave solution 

The basic properties of the system ([5]) can be found from the stability analysis of 
the fixed points on ±oo. On the right hand side, as £^ —^ oo, the fixed point 
coordinates (or asymptotic values of u and v) are given as u = and v = 1. 
We will refer to this stationary point as the end point. The linearization of ([5]) 
around these values suggests that the end point is a saddle-node for all physical 
parameter values with corresponding eigenvalues given as fii — 0, 112 = — c, 



In the opposite limit ^ ~oo, the fixed point (referred to as the starting 
point) has coordinates given by ([5]), i.e. any point lying on a line v = 1 — (3u 
in (u, v) plane can be a fixed point (or correspondingly, asymptotic values of 
temperature u and fuel concentration v belong to this line in the (u, v) plane as 



u = 13 




V = a, 



(6) 



fi3,4 = (-c± v/c2 + 4r/3)/2. 
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<^ — + —oo). On the (/3m, v) plane a travelling wave solution can be schematically 
represented as shown in figure[TJ The linearization of the governing equations ([5]) 
around the values (O suggests that the starting point can be either a saddle-node 
or a stable node depending on the parameter values. Namely, the eigenvalues 
of the linearized problem can be written in this case as: fii = 0, /Lt2 = — c, 
= i—c± \/c? + 4/3(r — fTe~''/(i~'^)))/2. The starting point is a saddle-node 
for T > cre"^/'^"'^"''^-' and a stable node for r < ae~^/'^^~'^\ In the latter case it 
is impossible to connect the starting point (as ^ — > —00) with the end point (as 
^ 00) with a trajectory in the phase space and therefore the solution does 
not exist. The condition 

r = a exp (7) 

1 — (T 

defines a surface in the (r, P, a) parameter space which represents the boundary 
between the region, where the solution exists from the region where it does not. 
If any of the solution branches crosses this surface it has to exhibit extinction. 
Therefore, we can also refer (O as the extinction condition. In figure 12.21 the 
boundary between the parameter regions where the travelling solutions exist 
(above the surface) or cease to exist (below the surface) is shown. The physical 
meaning of this condition can also be explained if we return to equation ([1]) 
and neglect the diffusion terms. Behind the travelling wave there has to be a 
stable state defined by ^ and w = 0, where no radicals can be produced if 
some small variation of radical concentration is imposed. As seen from ([T}, the 
latter condition is satisfied if r > a exp and the right hand side of the last 
equation in ([1]) is not positive. 

It is useful to consider several cross sections of ([7]) with the planes r = 
constant, since we will be studying the properties of the solutions to ([5]) for 
fixed values of r. The corresponding level curves are shown in figure [3] and are 
given as /3e(c) — (1 — cr) ln(cr/r). It is seen that /3 as a function of a has a single 
maximum, /3m, for a fixed value of r. We will denote the value of the residual 
amount of fuel for which /3(ct) has maximum as Cm, i.e. /?„ — (3{am)- This 
maximum is a solution to the d(3e/da — equation which can be written as 
(1 — CTm)/o'm — ^Tn.am/r — using equation ([7]). The solution of this equation 
is shown in figure S] and corresponding value /3m as a function of r is plotted 
in figure [5] It worthwhile noting that according to [15] the condition for the 
extinction can be expressed using our variables as /3 — In r and it is shown 
in figure [H] with a dashed line. The asymptotic expression for the extinction 
condition was derived in the limit /3 — > 00, when a 0. For moderate values of 
the activation energy this dependence is given by ([7]) with the slope 1/(1 — a) 
times smaller. Therefore the solid curve representing our results lies below 
the dashed line corresponding to the asymptotic prediction of [15j . For r — » 1 
parameter am tends to 1 and the curve /?,„ vs. logr becomes tangent to abscissa 
axis. In the opposite limit r ~* parameter /3,„ tends to 00 and cr„i to zero, so 
that the slope of both curves becomes equal. 

It is interesting to compare the burning temperature which we define as 
Uf, = (1 — ct) / f3 with the so called crossover temperature, Uc, a temperature at 
which the rate of branching and recombination are equal given that v — 1 i.e. Uc'. 
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Figure 2: The extinction condition ([7]) plotted as a surface in the parameter 
space (r, (3, a). 

r = exp(— 1/wc). Using this equation and relation r = cr,„ exp (—/?„/ (1 — <Jm)) 
we derive = Pm/i,^ — cfm) — lii(crm)- This yields 

Ub _ l-(T f Pm ^ \ _ Prn 1 - ^ 1 - (1 + In (Tm) ^-g^ 

Uc P \l - am ™/ /3 1 - CTm 1 - Cm 

Since P < Pm, a < Cm and In dm < the first two factors in the right hand side 
of the above equation are greater or equal to one and the last factor is always 
greater than one. This indicates that the actual temperature in the combus- 
tion wave is always above the crossover temperature (even for the extinction 
condition) . As discussed in [T3] the temperature of the flame has to be higher 
than the crossover temperature so that the rate of the radicals production can 
balance both their consumption due to chemical reaction and to depletion of 
the radicals due to diffusion. 

3 Travelling front solution 

In our previous paper [201 we have already shown that equations ([5]) subject 
to the boundary condition fl]) exhibit travelling front solutions similar to those 
found for the one-step reaction models, such that m(^) and w(^) are monotonic 
and and are bell-shaped functions of the spatial coordinate. Typical 
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Figure 3: The extinction condition ([7]) plotted as level curves for r — 0.1, 0.01 
and 0.001. 
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Figure 4: Dependence of maximal value of the residual amount of fuel, a„n on 
r for the extinction condition ([7|). 



10 



20 




log r 



Figure 5: Dependence of Pm on r for the extinction condition. The sohd hne 
represents the data plotted for ((7|) , whereas the dashed hne represents the results 
obtained using the paper [TS]. 

solution profiles w(^) and w(^) are shown in figures [S] and [T] The scaled coor- 
dinate X is used instead of ^ so that x S [0, 1]. As seen from the figures the 
solution profiles are sharper for smaller values of f3 and flatter for greater values 
of (3. The residual amount of fuel a for fixed r increases when /3 increased. 

The properties of the classical travelling front solution are summarized in 
figures [5] and [51 where the speed of the front c and the residual amount of fuel 
(T are shown as functions of /3 for several values of r. 

As noted in [20], at first glance the dependence of c upon /3 in figure [8] resem- 
bles the behaviour of the speed of the front for the model with one-step reaction 
scheme, which was studied in Namely, c reaches the global maximum for 
some value of f3 and decays monotonically as we increase (or decrease) /3 from 
the value corresponding to the maximum. However more detailed investiga- 
tion shows that there is a significant difference between the prediction of the 
one- and two-step models. In particular, for the model with one-step reaction 
mechanism the travelling front solution exists for any value of /3 and decays 
exponentially to zero as we increase /3. This is not the case for the model with 
two-step reaction mechanism studied here. As we increase /3 (for a fixed r) up 
to some critical value /3c, the speed of the front decays rapidly and the travelling 
front solution ceases to exist for [3 > (3c- Furthermore there is some residual 
amount of fuel left behind the front in the case of the two-step model unlike 
the one-step adiabatic model which has zero residual amount of fuel. To some 
extent, the properties of the two-step adiabatic model studied here resemble the 
properties of the nonadiabatic one-step model investigated in jl8J, more than 
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X 

Figure G: Classical travelling front solution profiles m(^) and v{£) for /3 = 0.054, 
r = 0.01. Solid line shows the temperature profile and the dashed line represents 
the fuel concentration. 
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Figure 7: Classical travelling front solution profiles u{S) and v{£) for (3 = 2.364, 
r = 0.01. Solid line shows the temperature profile and the dashed line represents 
the fuel concentration. 
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Figure 8: Speed of the front as a frinction of /? for various values of r. 
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Figure 9: Dependence oi a on (3 for various values of r. 
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the adiabatic one-step studied in [1]. This is expected since the recombination 
step, which is an inhibitor of the chain branching reaction, plays a similar role 
as the heat exchange with the surrounding in the one-step nonadiabatic model. 
In both cases there is a nonzero residual amount of fuel left behind the reaction 
zone. The similarity between these two cases is also strengthened by the like- 
ness of the behaviour of the speed of the front as a function of the parameter /3. 
Namely, in both the one-step nonadiabatic and the two-step adiabatic models, 
the travelling front solution ceases to exist as we approach some critical value 
of f3c (in the combustion literature this event is usually called extinction [TB]). 
However, the route to extinction in these models appears to be different. In the 
case of the one-step nonadiabatic model, for given parameter values, there are 
either two solution branches with different speeds or no solutions. The extinc- 
tion occurs when the two solution branches coalesce (this event is also known 
as a turning point or a fold bifurcation). For the two-step reaction mechanism 
the extinction occurs when the speed of the front drops down to zero as we 
approach the critical parameter values. 

It is evident that the route to extinction distinguishes this model from its 
well known one-step counterparts. In the next section we shall investigate the 
important phenomenon of extinction in greater detail. 

3.1 Route to extinction 

Previously we mentioned that extinction may occur when the parameter values 
of the problem reach the boundary of the existence of the solution ([7]). Does 
this imply that the the extinction of the classical travelling front solution occurs 
when the parameter values of the problem approaches the surface ([7])? Extensive 
numerical analysis shows that for fixed r the travelling front solution exhibits 
extinction as a and /3 reach the values (Jm and (the parameter values where 
the curve (3e{o'), which is an intersection of surface ([7]) with plane r = const ^ 
folds or the maximum value of (3 is reached for fixed r) . This is illustrated in 
figure [To] where the dependence oi a on f3 is plotted for r — 0.01 with a thick 
solid curve. The critical parameter values for the extinction are also shown on 
the same graph with a dotted curve. As seen the travelling front solution ceases 
to exist as f3^ Prn- 2.4153239 and a ^ 0.23947057. 

Next we investigate the behaviour of the travelling front solution for param- 
eter values close to the extinction point. Firstly, we note that as we approach 
the extinction point along the travelling solution branch, /3 [3m, c Cm, 
c ^ 0, ae~^/'^^~'^'^ r. In order to investigate the problem further we split the 
solution profile into several zones according to the physical processes dominat- 
ing in these regions. The region in front of the travelling wave, which is usually 
referred to as a preheat zone, is governed by two processes: diffusion and convec- 
tion. For this reason it is also sometimes called the diffusion-convection region. 
In this region the reaction terms are negligible and can be omitted from the 
governing equations. Without loss of generality let us assume that the preheat 
region satisfies f > i^f,. Since the travelling wave solution has a translational 
invariance, ^f, can be chosen to be arbitrary. For ^ > the front satisfies the 
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Figure 10: Parameter values a vs. /3 for the travelling front solution with 
r — 0.01. Thick solid curve shows the numerical results. Dotted curve represents 
the extinction curve i.e. intersection of surface ([7]) with plane r — 0.01. Dashed 
line marks the value am at which the extinction curve reaches the maximum 
value of p. 
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following equation 



+ CM^ = 0, 

ujj + cw^ — 0, 



(9) 



and boundary conditions 

M = 0, v^l, for ^ ^ 00. (10) 
A solution to ^ subject to ([T0|) can be written explicitly as 

u{O^Ube'<, v{0 = l-{l-Vb)e--^, (11) 

where Ub and have to be found by matching the solution with those of ^ < ^b- 
In the region ^ < ^6 the governing physical processes which occur in this zone are 
reaction and diffusion, and therefore is sometimes called the reaction-diffusion 
zone. The convection process is negligible when compared to other terms and 
can be omitted from the governing equations. Re-scaling the temperature as 
u — > (3u, spatial coordinate as z — > ^/firS, and dropping the convection terms we 
can rewrite the governing equations as 



(12) 



(1-u-i;) =0, 

'"-{l-u-v)e-P/'' 
r 

which has to be solved subject to the boundary conditions 

u=l — (7, V = a, z ^ —00. (13) 

In the opposite limit z — > zt we require that the reaction terms vanish as z — > z^, 
which implies immediately that u{zb) + v{zb) — 1 (or w = I — u ~ v ^ 0). 
Moreover, since the reaction terms are neglected in the preheat region we assume 
that derivatives of u{z) + v{z) also vanish at z = Zb. The parameter Zb is 
undefined and we shall take it to be a sufficiently large number and replace 
the problem on the half line z G [— 00, Zb] by the problem on an infinite line 
z G [—00, 00] which is more convenient for further investigation, noting that 
this expansion of the domain does not impact on our results. 

It is convenient to define new variables w = 1 -I- aui and v = avi and use 
the substitution r = a,^e~^"^/^^~'^"^^ to rewrite the governing equations as 



Ulzz ~ Ul - Wl = 0, 

Vlzz + — ^;i(mi -I- Vi) exp[/3m/(l - am) - /3/(l - (TUi)] = 0, 
CTm. 



(14) 



and boundary conditions as 



ui = —1, vi = 1, for z — !■ —00, 

ui+ vi = 1, uiz + Viz — 0, for z — s- 00. 



(15) 
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The solution to the problem P^HTS)) is sought in the form of the series 



Ui 



(z) = -l + ^^'4')(z), i;i(z) = l + X]5^i;«(z), (16) 



with S being a small parameter of the asymptotic expansion. The physical 
interpretation of this parameter will be defined later. As a next step we expand 
the exponential temperature dependent factor in the second equation of (jl4|) in 
a Taylor series about ui = — 1: 

(17) 

where the notation A ~ (t/3/(1 — cr)^ has been introduced and terms higher than 
quadratic in ui + 1 ^ 0{6) have been neglected. Substituting this expansion 



into IHl) yields 

Ulzz - ui - Vl = 0, 

vizz + Kvi{u, + V,) (l + A{u, + 1)+a(a- 2-^] (!fl±lL ) ^ q, 



(18) 

where parameter K — a /dm exp[/3m/ (1 — CTm) — /3/(l ~ f)] has been introduced. 

Since it is assumed that the parameters are chosen to be close to the ex- 
tinction values, we can therefore write a = am + Ac, where Ac 1. At 
the same time we know that /3 as a function of a is tangent to the extinc- 
tion curve I3e{d) at a — Um- Therefore it is reasonable to conclude that 
j3 « (3m + /52(A(t)^ -1- ... since (3e{<j) has a maximum at cr = Om- The coef- 
ficient (32 is equal to 1/2{(P (3 /da'^)\a We shall now proceed to show the 
dependence of K and A on Act. By differentiating ([7]) with respect to a for fixed 
r it can be shown that cr/?e((T)/(l — cr)^ ^ 1 as ct — > am- Therefore, it can be 
written that A « 1 + AiAtr + 0((Acr)2), where Ai = /?,„(! -I- tTm)/(l - amf ■ 
On the other hand using the expansion for (3{a) and A it can be shown that 
K^l- K2(Acr)2, where = (1 + ct™ + 2(32al,) / {2al^{l + amf)- It is conve- 

1/2 

nient to use the small parameter defined by the equation 5 = Acrifj instead 
of Act. In this case we can write the following asymptotic expansion for the 
parameter values near the point of extinction 



= CT„i + A2 0, 
(3 = Pm+(32K^^5^ + - 

A = 1 + Ai5 + 



(19) 



where Ai — A^jK^ 

Next we reduce the order of the differential equations (fT8| by redefining the 
new independent variable as v instead of z (here and later on, it is convenient 
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to drop the subscript '1' in (HI])). This implies immediately that we restrict 
the consideration only to the case of solutions with monotonic dependence of v 
on z. However, our travelling front solution always satisfies this property, and 
therefore we do not impose any restriction by replacing z with v. Introducing 
p = Vz, and substituting d/dz = pd/dv and d'^/dz'^ = pp^d/dv +p'^d^ /dv"^ into 
pg]) yields 



P^Uyy + pPvUy - U- V = 0, 



(7 \{U+1) 



(20) 



ppy + Kv{u + v) \ l + A(m + 1) + A - ^ j ^ — 7^ ) = 0, 

subject to the boundary conditions 

u = —1, P = 0, for u 1, 

u + v — Q, Uv = —\, for u — > oo. 



(21) 



Earlier in this section we have noted that the travelling front solution is sought 
in the form of an asymptotic series (see equation [T6|) so we can expect that the 
solution profiles change substantially on the length scale of the order of 5 in v. 
On the other hand, it can be shown that close to w — 1, the solution to (I20p 
behaves as u ^ — 1 — K{v — 1), u + w ~ 5'^{v — 1), and p ~ 6{v — 1) i.e. u, u + v 
and p are of different order in terms of the small parameter 5. In order to take 
this into account we introduce a new variable s{v) = u + v and seek the solution 
in the form of 

s(v) = 5^S3,{v) + 5^Si{v) + S^s^iv) + , 
p{v) = S^P2{v) + S^Psiv) + 6'^Pi{v) + ... ^ ' 

We also re-scale the independent variable in ((20|) asv = l+(5y, so that y e [0, oo). 
Substituting (|22p and the parameter expansion (fT9|) into (|20|) and using ?/ as a 
new independent variable we obtain a system of differential equations to each 
order of the small parameter 5. Combining these equations so as to leave only 
the leading order terms S3 and p2 yields a closed system of equations 



(23) 



phi + P2P'2S3 - sz{l + \iy + ny^) =0, 
P2P2 + S3 = 0, 

subject to the boundary conditions 

P2 = 0, S3 = 0, for y ^ 0, ^24) 

P2 = const, S3 = 0, for y 00. ^ ' 

Here we have introduced k = (1 + (Tm)/2(1 — Um) and the prime stands for 
d/dy. The leading order problem ([23l - [24|) includes two parameters Ai and n. 
The latter is a fixed given number for each value of r, whereas the former is an 
eigenvalue which has to be found together with the unknown functions P2{y) 
and S3(7;). The constant value V is defined from the solution P2{y)- 
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Once the problem (|^5H^ has been solved for a fixed r, the parameter values 
Ai, K2 and (32 can be found from their definitions above. The only unknown pa- 
rameter is the speed of the front, which can be found by matching the solutions 
in the reaction-diffusion and preheat zones. In order to undertake this matching 
we have to transform the variables P2(,y) and s^(y) back to the original vari- 
ables u{£) and v{£) which are used in equations (O. To summarize the various 
sequences of variable changes we provide the following list of transformations 

Jm (3(1 _ a{s{y) - 1 - <5y)), ^35) 
^El ^Gl crvi(z)) Gv a{\ + by). 

The index above the arrows shows the equation number, where the correspond- 
ing substitution has been made. In the leading order 0((5°) we have the match- 
ing condition: u\y — \ — a and vy, = cr, where u\y and Vh are taken from 
We also have to match the derivatives of the solutions to the boundaries of the 
reaction-diffusion and preheat regions. Solution of the preheat problem suggests 
that — c?w/d^|b — dv/d^\}y = (1 — cr)c. On the other hand, it follows from the 
boundary conditions ([M)) that 



dv 



b 



= cr,n\/ Pmr hm p = 

(26) 

Therefore, we have the following expression for the speed of the front 



c = ^M^a.mVSy{l - a„0, (27) 

where terms of the order higher then 0(6^) have been dropped. 

At a first glance, the asymptotic procedure of the speed derivation is incon- 
sistent. Indeed, the convection term in equations ()12|) which is proportional to 
c ~ 5^ has been neglected. Nevertheless, from equations ([T51) onwards the terms 
of the order of 6^ are retained. This, however, does not lead to any inconsisten- 
cies because if the convection terms are retained throughout the analysis in the 
reaction-diffusion region, they still will not effect the leading order equation (|23p 
and will appear only in the equations for the higher order of small parameter 
6 (i.e. in the equations for p^ and 54 etc.). Therefore in terms of the original 
parameters of the problem we have 

/3 = /3m(r)-/32(r)(a-a™(r))2, 

c = ^/^-^^V{r)K2ir){a - a,„(r))2. (28) 
t — Cm 

In order to obtain estimates for /3(ct) and c(ct) using equations (j28p we are 
required to solve the system ((23)) numerically to obtain Ai and V. The solution 
profiles P2{y) and 53(2/) are shown in figures fTTl and fT2l For small values of y, 
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2.0 




Figure 11: Solution profiles P2{y) for equations ([^5)1 for various r values. 



both p2 and S3 tends to zero linearly as y and —y respectively. In the opposite 
limit of large y the value of S3 tends to zero. This causes the reaction terms 
in ([221) to vanish and p2 to reach a constant value which defines the speed of 
the front and the derivatives of the temperature and the fuel concentration at 
the boundary between the preheat and reaction-diffusion regions as discussed 
above. Once the system ([23]) is solved for some r, the dependence of the front 
speed c and /? on the residual amount of fuel, cr, can be found using ([28]) . In 
figures [12] and [TH c and fi are shown as functions of Act — d — tim for the 
values of a close to the extinction value am and r = 0.1. The numerical results 
obtained by solving ([5]) are presented on both figures with a solid line, whereas 
the dotted lines shows the prediction of the asymptotic analysis (|28|) . In figure 
[13] we also plotted the extinction curve /3e(c) with a dashed line. As seen from 
these figures the correspondence between the numerical and asymptotic results 
is excellent for a am. As |Act| increases up to the value of 0.1, the higher 
order terms in the asymptotic expansion begin to play a significant role and the 
discrepancy between the numerical and asymptotic results becomes visible. 

We also compare the numerical and asymptotic results for different values of 
r. In figure [T5I and [TBI we plot the dependence of P2 and C2, which is 2(Pc/da'^ 
as a ^ am, on r respectively. Crosses connected with dashed line show the re- 
sults obtained analytically by employing (j28p and squares denote the numerical 
solution of ([5]) by calculating 2(P(3/da'^ as a ^ am in figure [15] and 2(Pc/da'^ 
as CT — > am in figure 1161 The correspondence between the numerical and ana- 
lytical results is excellent for r ranging over several orders of magnitude. The 
discrepancy found in both cases was only in the third significant digit. 

To conclude this subsection we shall summarize the main results obtained 
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Figure 13: Dependence of /? on ct for r = 0.1. Solid curve is plotted using the 
results of numerical integration of ([5]). Dotted curve represents the prediction 
of asymptotic formula ([28]) . Dashed curve shows the extinction curve (3e{<^) for 
r = 0.1. 
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Figure 14: Dependence of c on cr for r = 0.1. Solid curve is plotted using the 
results of numerical integration of Dotted curve represents the prediction 
of asymptotic formula 
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Figure 15: Dependence of /32 on r. Crosses connected with dashed line shows the 
results obtained with asymptotic formula ([28l) . Squares represent P2 estimated 
based on the numerical solution of ([S]). 
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Figure 16: Dependence of C2 on r. Crosses connected with dashed hne shows the 
results obtained with asymptotic formula ([28|) . Squares represent /32 estimated 
based on numerical solution of ©■ 
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Figure 17: Dependence of dc/d(3 on r. Triangles shows the results obtained with 
asymptotic formula (|28|) . Solid line represents the asymptotic results derived in 
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here. We have investigated the behaviour of the traveUing front solution for pa- 
rameter values close to the extinction point. We determined that the extinction 
of the travelling front solutions occurs when (3 and a reach the critical values 
Pm and am corresponding to the maximum condition for the function (3e{o') de- 
termined by ([7]) for a fixed r. The analysis of the solution behaviour around the 
critical conditions and <Tm shows that the solution branch on /3 vs ct follows 
a quadratic law which is also the case for the dependence of the speed of the 
front on a (see equation (^5]) ). The dependence of speed of the front vs. /3 is a 
linear function for /3 close to (3rn ■ Here we note that the same type of behaviour 
for c(/3) is predicted in [TSj , although in [IS] a is not considered and the authors 
employed different nondimensionalization. In terms of the parameters used in 
this paper the result of [15] can be written as c ~ y^r/ \n{l/r){Pm — /?) . In figure 
[TTI the derivative of c vs. (3 ai (3 = (3m obtained by using this formula (solid 
line) and calculated with asymptotic formula ([28]) (triangles) are compared. It 
is seen these results are in good agreement. In the next subsection we proceed 
to investigate the stability of the travelling front solution. 



3.2 Stability of the travelling front solution 

We investigate the stability of the travelling front solution numerically using 
two methods: direct integration of the governing PDEs and the Evans function 
method '21]. For both of these methods it is convenient to rewrite ([T]) using 
the assumption of equal diffusivity of fuel, heat and radicals (i.e. ta = tb = 1) 
and the integral relation S = f3u + v + w = I (which is employed to reduce the 
governing equations ([3]) to ([5])) as 

ut = Uxx + r{l- Pu-v), , . 

vt=vxx- Pv{l- f3u-v)e"^/^, ^ ' 

subject to the boundary conditions 

u — 0, V — 1 for a; ^ oo, , , 

Ux — 0, Vx — for X — oo. ^ ' 

In order to investigate the linear stability of the travelling front solution u{^), 
v{^) using the Evans function method, equations (j29p are linearized around 
it and we seek the solution of the form: u{x,t) = u(^) -I- £iy9(^, t)e'''* and 
v{x,t) = v{£^) -|- ex(^)e^*. Here (p{£,) and xiO ''^^^ linear perturbation terms, 
e is a small number, and A is the spectral parameter which is a complex num- 
ber. Substituting this solution into the governing equations (|29p and retaining 
terms linear in e yields 

^:)-^(:)^ (31) 



X J \ X 



where 



9| + — r/3 



^ ^(/3u2+/3u + v-l)e-i/" d1 + cd^ + (3{2v + I3u - l)e-^/'' 



(32) 
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The stability of the travelhng front is then defined from the spectrum of L. We 
can show that the essential spectrum of this operator always lies in the left half 
plane by using the approach outlined in [52] . We consider the limiting operators 
= limj_,±oo L. It can be shown that the essential spectrum of consists 
of algebraic curves 

Ai(fc) = +icfc, A2(fc) = -fc2-/3(r-(Te-'3/(i-'^))+icA:, , . 

where k e (—00, +00), which are located in the left half plane and are symmetric 
about the real axis. The first curve in ([551) includes the origin. Suppose that 
K is the union of the regions inside or on these curves, so that C\K contains 
the right half plane. Then according to [22], the essential spectrum of L is 
contained in K, and in particular includes curves ([33]) . Therefore the essential 
spectrum of L lies in the left half plane (including the origin) and the discrete 
spectrum is solely responsible for the transition to instability. Namely, it can 
be said that the travelling front is linearly unstable if, for some fixed complex 
A with i?e(A) > 0, there exists a solution of ((3T]) which decays exponentially 
as ^ — > ±00. We will refer to this A as an eigenvalue and to the corresponding 
solution as an eigenmode. 

We investigate the discrete spectrum of L using the Evans function method. 
The Evans function D{X) is a function of the spectral parameter A analytic in 
C\K. It has the following properties: D{X) = if, and only if, A is an eigenvalue; 
the order of any zero of the Evans function D{X) and the algebraic multiplicity 
of the corresponding eigenvalue are equal. Therefore, the stability analysis of 
the travelling front of ([29| can be reduced to the search for zeros of the Evans 
function located in the right half plane. In our previous papers [U [T5] we have 
shown how to construct the Evans function and how to calculate it numerically. 
In [3] it was demonstrated that the number of zeroes N located in the right 
half plane can be calculated using the Nyquist plot technique. Obviously, if 
N = then the travelling front is stable and if TV > then the travelling front 
is unstable. We have calculated N along the travelling front solution branch of 
(|29]) for various values of r = 0.1, 0.01, and 0.001. Parameter f3 has been taken 
for a range from /3 < 0.01 up to for each value of r. In all cases we have 
found that iV = 0, that is the travelling front solution is always stable. The 
same result was obtained in 15J in the limit of high activation energy. 

We also investigate the stability of the travelling front solution by direct 
integration of the governing PDEs ([55)1 using the time and space adaptive finite 
element package f23^ . We take the initial solution profiles in the form consistent 
with the boundary conditions (j30p as 

uix,0) = (2/3)-i(l - a)(l - tanh[lix - xq)]), , 
v{x,0) = (1 + cr + (1 -cr)tanh[;(x-a;o)])/2, ^ ^ 

where I and xq are chosen in such a way so that u{x,0) and v{x,0) fit the 
travelling front profiles u(^) and w(^) of ([5]). In figures [TSl and fT9l tvpical solution 
profiles u(x,t) and v{x,t) of ([5]) are shown. The initial conditions converge 
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Figure 18: Temperature profile u{x, t) obtained by numerical integration of (|29p 
for /3 = 1 and r = 0.01. 

quickly to a travelling front solution which propagate with constant speed. Note 
that the time axis in figure [19] is inverted in order to give a proper perspective. 

Solutions to (f29|) for different values of (3 and r = 0.01 are also calculated. 
The travelling front is found to be stable for all trial parameter values. In 
figure we plot the speed of the travelling front as function of [3 for r — 0.01. 
The numerical results obtained from direct integration of the system of ODEs 
([5]) (solid line) and PDEs (squares) are represented on the same graph. The 
correspondence between the prediction of both approaches is excellent. The 
discrepancy was found only in the fourth significant digit. 

It is interesting to study the behaviour of the solution of the governing PDE 
(PS]) for various values of (i above f3m in order to investigate if the flame can 
exist beyond the critical extinction condition. We integrate equations (|29p for 
r — 0.01 and (3 — 2.5 > (3m- For all trial initial conditions of the form discussed 
above, the solution decays to the equilibrium state u{x,t) = and v{x,t) = 
1 with typical solution profiles shown in figures [2T] and [22l The maximum 
temperature u{0,t) drops rapidly in time to zero and subsequently the reaction 
ceases and the reaction terms vanish in (|29)) and the system exhibits almost 
linear diffusion to the equilibrium state. This is clearly seen in figure E51 which 
shows the dependence of the squared inverse maximum temperature upon time. 
The numerical results obtained with the integration of the governing equations 
(|^^ (dashed line) are compared with the analytical prediction (solid line) , which 
was found by considering equations ([29]) without any reaction terms (i.e. we 
consider only the diffusion effects) with initial conditions as described above. In 
this case it can be shown that in the limit t oo the asymptotic behavior of 
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Figure 19: Fuel concentration profile v{x,t) obtained by numerical integration 
of ([291) for /3 = 1 and r = 0.01. 
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Figure 20: Speed of the travelling front solution as a function of (3 for r — 0.01. 
The solid curve represents the results of integrating the system of ODEs ([5]), 
whereas the squares correspond to the results obtained from integrating the 
governing PDEs ((29)) integration. 
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Figure 21: Temperature profile u{x, t) obtained by numerical integration of (|29p 
for (3 = 2.5 and r = 0.01. 

the maximum temperature satisfies u(0,t)~^ = ■nP'^t/ xf^{l — a)"^ , which is shown 
with a solid hne on the figure 1231 The numerical results are shifted from the 
asymptotic prediction, especially near the starting point < = 0. However for 
larger values of t the maximum temperature follows the diffusion law derived 
from the linear diffusion equations. 

4 Conclusions and discussions 

In this paper we have studied combustion wave propagation in an adiabatic 
model with two-step chain branching reaction mechanism. The current work 
extends our previous preliminary investigation of chain branching flames j20| . 
In contrast to [15j we are not using the activation energy asymptotic approach 
and consider the model for general values of activation energy. We have also 
used different non dimensional parameters which enabled a more convenient 
comparison between the properties of one- and two-step models. 

It is demonstrated that the model exhibits travelling combustion front so- 
lutions. The properties of these solutions differ form the properties of one-step 
models. Combustion waves are found to exist in certain regions of the param- 
eter space and are characterized by non zero residual amount of reactant, a, 
left behind the travelling wave. This is not possible in one-step adiabatic flame 
models as all the fuel is used up in such models. To some extend this charac- 
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Figure 22: Fuel concentration profile v{x, t) obtained by numerical integration 
of ([29]) for /3 = 2.5 and r = 0.01. 
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Figure 23: Decay of the maximum temperature u(0, t) of (|29p solution with time 
for r = 0.01 and (3 = 2.5. Dashed line represent numerical results and solid line 
shows the analytical estimate of decay rate. 
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teristic of the two-step adiabatic model is similar to the properties of premixed 
combustion waves in nonadiabatic one-step models, which can also exist in a cer- 
tain region of parameter space and exhibit extinction as the boundaries of this 
region are reached. However the behaviour of the travelling combustion waves 
near the extinction condition are completely different in these two types of mod- 
els. In the one step nonadiabatic model the point of extinction corresponds to 
the turning point of the dependence of flame velocity on parameters (which is 
a two- valued function), therefore the speed of the wave is always positive. In 
our two-step model studied here the flame velocity is a single valued function 
and the velocity drops down to zero as the extinction condition is reached. It 
is remarkable that at the extinction condition there exist a flame characterized 
by stationary distribution of temperature, reactant and radicals i.e. a standing 
combustion wave. 

The boundaries of existence for combustion front are determined. We show 
that for a fixed value of r (dimensionless recombination rate) travelling wave 
solution branch occupy a certain region in the parameter space which lies outside 
the region bounded by the curve /9e(c) which has a maximum /3m at am so that 
P remains less than or equal to Pm and a residual amount of fuel less than am- 
The maximum values f3m and achieved at the point of extinction. 

The properties of the travelling front solution are studied both analytically 
using the matched asymptotic analysis and numerically via the integration of the 
governing PDEs and the ODE formulation of the travelling wave solution 
(O. The matched asymptotic analysis allowed us to obtain the properties of the 
travelling front solution near the point of extinction. The asymptotic expression 
for the speed of the front is obtained as well as the location of the solution branch 
in the parameter space. In particular, it is shown that for a fixed value of r the 
speed of the front and f3 are quadratic functions of a — am near the point of 
extinction. The results of the asymptotic analysis are then compared with the 
results obtained by the integration of the governing ODEs for the travelling wave 
solutions by the relaxation method as described in [H [TS] . It is demonstrated 
that the results obtained from these two approaches are in excellent agreement. 
The speed of the front and the location of solution branch in the parameter space 
are studied in detail for a broad range of parameter values by means of numerical 
integration of the ODEs. The results of this investigation are then compared 
with the direct numerical integration of the governing PDEs and are shown to 
be in an excellent agreement. The stability of the travelling front solution is 
investigated numerically by employing the Evans function method. It is shown 
that the classical travelling front solution is stable for all parameter values. 
These results are also confirmed by direct numerical integration of the governing 
PDEs. We also compare the results of our paper with the results obtained in 
|15| in the limit of high activation energy and show that they qualitatively agree. 

In this paper we considered the case when the diffusivity of fuel and radicals 
are equal: t — ta = tb- Moreover, for the sake of simplicity it is assumed 
that these diffusivity are both equal to the heat diffusivity: r = 1. Using 
these assumptions we found that the classical travelling solution is stable for 
the whole range of parameters r, /3, and a where the solution exists. Similar 
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situation was observed for the one-step nonadiabatic flames. In [T5] it was 
shown that for r = 1 one of the traveUing front solution branches was stable 
for all parameter values where the solutions exist. However, for different values 
of r < 1 the situation changes due to the presence of the Bogdanov-Takens 
bifurcation. Therefore, one of the most important issues for further investigation 
is to understand how the stability of the solutions changes in the chain-branching 
model for different values of t. 
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